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Abstract: We investigate two common approaches to model development for robust control 
synthesis in the aerospace community; namely, reduced order aeroservoelastic modelling based 
on structural finite-element and computational fluid dynamics based aerodynamic models and 
a data-driven system identification procedure. It is shown via analysis of experimental Super- 
Sonic SemiSpan Transport (S 4 T) wind-tunnel data using a system identification approach it is 
possible to estimate a model at a fixed Mach, which is parsimonious and robust across varying 
dynamic pressures. 

1 INTRODUCTION 

Undesirable aeroservoelastic (ASE) interactions are a major concern in modern aircraft design. 
ASE interactions between aircraft structure, aerodynamics and flight control systems can lead 
to divergent oscillations resulting in catastrophic failure [1], As such, analytical model devel- 
opment is an important step in the design and certification of aircraft. Accurate models allow 
for robust control design, which is critical for aircraft safety, gust-load alleviation, ride quality, 
etc. 

Finite element (FE) based models are used in the design process of aircraft to aid in the de- 
scription of complex elastic and structural components [2,3]. FE based models that accurately 
characterise the aerodynamic and structural components are of very high order (e.g. thousands 
of degrees of freedom) and computationally intensive. To reduce the ASE model order a modal 
approach is used. This can reduce the state-space model order to several tens of states. The 
robustness of this model type highly depends on the accuracy of the FE structural and aerody- 
namic models and on the number of states applied to the modelling. Models of high complexity 
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inhibit their use for control synthesis because their real-time implementation is difficult or not 
possible [4]. This difficulty has led to considerable activity in the areas of model and controller 
reduction techniques in the last decade. The literature is rich with many reduced order model 
(ROM) techniques which would require a review paper to properly discuss them. As such, in 
this paper, we limit ourselves to a handful of approaches. 

A leading strategy to FE based model reduction is the proper orthogonal decomposition (POD), 
which is also known as the Karhunen-Loeve procedure [5]. The so-called POD technique is 
well known in the statistical literature as principal-component analysis [6]. The reduced basis 
method was first proposed in [7, 8] for structural analysis and it has been used for structural 
problems in [9-11]. This technique uses synthetic data from a high fidelity FE based model 
to capture the dominant characteristic information utilising an orthogonalisation process. This 
allows the POD approach to accurately describe a system using a few basis terms, which gives 
it an advantage compared to other numerical procedures [12]. These reasons have allowed POD 
to become a popular technique for the implementation of real-time control [13]. Moreover, it 
has been successfully used in a variety of fields including signal analysis and pattern recognition 
[14], fluid dynamics and coherent structures [15-17], control theory [18-20], civil engineering 
[21] and inverse problems [22]. 

More recently the ROM techniques utalising POD were developed for aeroelastic systems anal- 
ysis [23]. This methodology was introduced to the aerospace community for the reduction 
of aeroelastic equations [24]. As a follow-on to this technique, frequency-domain approaches 
were developed which efficiently compute POD basis functions for linearised aeroelastic sys- 
tems [25, 26]. Due to its popularity and utility POD methodology has been proposed and im- 
plemented for static and dynamic continuous-time nonlinear aeroelastic problems. Subsequent 
development lead to extensions to encompass discrete Euler equations [27,28]. Successful ap- 
plication of this approach has been demonstrated for the analysis of limit-cycle oscillation of an 
airfoil with a nonlinear structural coupling in the transonic regime [29]. 

Nevertheless, it has been observed that standard POD procedures are less robust for nonlinear 
problems and typically require more basis functions as the function complexity increases [23, 
30]. To address this limitation a linearisation strategy, the trajectory piecewise linearisation 
(TPWL), was proposed [31,32]. The TPWL technique combines reduced-order modelling with 
linearization of the governing equations as a solution to this problem [23]. 
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Although POD offers a significant reduction of the full FE based model, it is often too large to 
lend itself to efficient control design. An alternative approach is to use data driven techniques 
to let the data dictate what the optimal model should be. Such a procedure is commonly known 
as system identification. 

This area, as with ROM methods, also has an extensive base of literature for both linear and 
nonlinear system identification techniques and is too large to give a proper review here. As 
such we refer the reader to an often cited authoritative treatise in the area which provides an 
excellent overview and references [33]. Below we provide a brief and incomplete introduction. 

There are two broad classes of techniques that can be pursued to accomplish the task of sys- 
tem identification: (i) nonparametric and (ii) parametric methods. The finite impulse response 
function (FIR) has been widely used for modelling linear time-invariant systems. This type 
of system description is known as nonparametric because it is a numeric representation of the 
system’s impulse response or kernel [33-35]. Although nonparametric methods can be used to 
represent many classes of systems, they do so at the expense of introducing an excessive num- 
ber of unknown coefficients which must be estimated. Most expansions map the past inputs 
into the present output and so require a very large number of coefficients to characterize the 
process. Moreover, the parameters are not readily linked to the underlying system, except in 
special cases where significant a priori knowledge of the system has been assumed. 

Due to this shortcoming, parametric identification methods have been developed for use in 
the design of better control systems. Parametric models have some advantages in applications. 
They (i) are easier to understand and interpret, (ii) can simplify forecasts (e.g., obtaining forecast 
intervals) and (iii) model comparison in a parametric context (i.e., parameter estimates, model 
order and model structure) has been well studied. Hence, the difficulty of model comparison 
encountered using nonparametric tools can be avoided [36]. 

A parametric model consists of a set of differential or difference equations describing the system 
dynamics. Such equations usually contain a “small” number of parameters which can be varied 
to alter the behavior of the equation. In this paper, we only consider the discrete-time case since 
in any practical experimental situation the data available to the experimenter is in discrete-time. 
As such, most systems for identification purposes are represented in discrete-time. In addition, 
we assume the ROM with which we compare our data-driven model is available and as such we 
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do not discuss its development. 


The organization of this paper is as follows. In §2 we formulate the identification problem 
addressed here. Section 3 describes the experimental S 4 T test-bed [37] and methods used for 
model development. Section 4 illustrates the results of our study on three flight conditions; 
namely, model robustness to fixed Mach but varying dynamic pressure. Section 5 provides a 
discussion of our findings and §6 summarises the conclusions of our study. 

2 PROBLEM STATEMENT 

System identification is the process of developing or improving a mathematical representation 
of a physical system based on observed data. Often the observed data consists of an external user 
selected input, used to perturb the system and elicit an output response. In any experimental 
situation the system output is a sum of the true unknown system output and observation or 
measurement noise. Given this paradigm there are several model structures that can be explored 
to model a system’s dynamics and noise. 

2.1 Model Structures 

After preprocessing the recorded data, the first step in system identification is to select a model 
structure to describe the observations. There are a number of model forms to select from when 
developing a data-driven model [33]. However, the problem often dictates the model form(s) 
that are reasonable to consider. This insight may come from a priori knowledge of the physical 
process or previous morphological modelling efforts. In this paper, we use knowledge gained 
from previous work to limit the model sets considered, namely, linear, time-invariant processes 
where both input-ouput are available to the user [38,39]. 

2.1.1 ARX Model Structure 

The simplest input-output polynomial model is the AutoRegressive exogenous input (ARX) 
model, represented as [33,40,41] 

y(n) = - aiy(n - 1) a na y(n - n a ) (1) 

+ b\u{n — 1) + b 2 u(n — 2) + ■ • • + b nb u(n — rib) 

+ e(n) 
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where y(n) is the measured output, u(n) is an accessible input, and e(n) is an unobservable 
white-noise disturbance. The current output depends on rib previous values of the input, n a 
previous values of the output and the current disturbance. 

This structure can be represented more compactly as 

A(q)y(n ) = B(q)u(n ) + e{n) or (2) 


y(n) = G(q)u(n) + H(q)e(n ) where 


(3) 


G(q) = 


B(q ) 


H(q) = 


(4) 


A(qY A(q) 

where A(q) = 1 + dig' 1 + ■ • • + a Ua q~ na , B(q) = b x q ^ + • • • + b nb q~ nb , g” 1 is the backward 
shift operator and the a’s and 6’s are the parameters of the output and input, respectively. 


When there is evidence of significant noise in the system a more flexible noise model may be 
required to model the dynamics and noise using different polynomials. 


2.1.2 ARMAX Model Structure 

For linear systems, the relationship between input-output and noise can be written as a linear 
difference equation 

y(n) = - cti y(n - 1) a Ha y(n - n a ) (5) 

+ b\u{n — 1) + b 2 u(n — 2) + • • • + b nb u(n — rib ) 

+ e(n) + cie(n — 1) + F c„ c e(n - n c ). 

This is known as the AutoRegressive, Moving Average exogenous (ARMAX) model. In this 
model structure the current output y(n) depends on an exogenous input, u(n), an innovation 
process, e(n), and past values of the output. This structure can be represented more compactly 
as 


A(q)y(n) = B(q)u(n) + C(q)e(n). 


Substituting (6) into (3) yields 


G(q) 


B(q ) 
A(qY 


H{q) 


C(q) 

A(q) 


(6) 


(7) 
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where C(q) = 1 + cig” 1 H h c nc q~ Uc , and the c’s are parameters of the noise model numer- 

ator. The extra polynomial, C(g), gives the ARMAX structure additional flexibility to model 
additive disturbance. When additional complexity is needed to model noise H (g) can be fully 
parameterised independent of the system dynamics. 


2.1.3 Box-Jenkins Model Structure 


A natural development of the ARMAX model structure is to parameterise the noise process as 
an ARMA model 


y(n) 


B(q) C(q ) 

M n ) + 7yrw e ( n ) 


F(q) 


D(q) 


(8) 


where F(q) = 1 + /ig” 1 + • • • + f nf q~ nf and D(q) = 1 + dig” 1 + • • • + d nd q~ nd model the 
poles of the system and noise separately. Note that rif — n a . This model structure is known as 
the Box -Jenkins model due to their seminal work proposing this model form [42]. 


2.2 Model Order Selection 

Once a model structure is selected the dynamic order of the system needs to be determined. 
Although there are many techniques that offer a solution to this problem, below we describe 
three commonly used techniques to estimate model order. 


2.2.1 Final Prediction Error 

The Final Prediction Error (FPE) measure estimates the error in model fit when it is used to 
predict new outputs [43]. FPE defines an optimal model as one that minimises 

where N is the number of data points, V is the prediction error, or the residual sum of squares 
and p is the number of model parameters. 

2.2.2 Akaike’s Information Criterion 

Akaike’s Information Criterion (AIC) is a weighted estimation error based on the unexplained 
variation of a given time series with a penalty term when exceeding the optimal number of 
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parameters to represent the system [44]. Utalisng AIC, an optimal model is defined as one that 


minimises 

AIC = v(l + ^y (10) 

According to Akaike’s theory, the most accurate model has the smallest prediction error. 

2.2.3 Minimal Description Length 

Rissanen’s Minimal Description Length (MDL) approach is based on V plus a penalty for the 
number of terms used [45]. With MDL, an optimal model is one that minimises 

( p In N \ 

MDL = Wl + C— ). (11) 

A model that minimises the MDL allows the shortest description of measured data. 


2.3 Cross-Validation 


Model validation is an important step in developing strategies for robust control. This step 
is typically preceded by system identification. Model validation is concerned with assessing 
whether a given nominal model can reproduce data from a plant, collected after some initial 
experiments to obtain estimation data [46]. The model validation problem is really one of 
model invalidation since a given model can only be said to be not invalidated with the current 
evidence. Future evidence may invalidate the model. 


We cross-validate the parameter estimates of the model dynamics using a 1-step-ahead predictor 
[33]. Model goodness is assessed by computing the percent quality of fit (%QF) as 


%QF = 


1 - 


Vn 


X 100 


where y is the measured system output and y the predicted model output. 


( 12 ) 


In the sequel, we show by implementing these well-known identification strategies it is possible 
to develop a model which is parsimonious and a robust predictor of measured data and, hence, 
represent the physical process more accurately. 
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Figure 1: S 4 T model. 

3 EXPERIMENTAL S 4 T WIND-TUNNEL DATA 

The modelling and identification techniques were applied to experimental wind-tunnel data 
from the S 4 T project conducted at NASA Langely Research Center [37]. Figure 1 shows a scale 
model of the S 4 T testbed in the wind tunnel. The data analysed for this study used horizontal 
tail position input and structural accelerometer response output. 

3.1 Data Collection 

Wind-tunnel data was gathered during transonic clearance of the S 4 T. At each flight condition 
the aircraft model was perturbed with a log sine sweep input which had a frequency content 
of 0.5-25 Hz, mean value of 3.5 deg and ±0.3 deg amplitude. The inputs were applied to the 
horizontal tail. The accelerometer output was collected from a sensor located at the nacelle 
inboard aft position. Wind-tunnel tests were conducted at subsonic, transonic and supersonic 
conditions and varying dynamic pressures. The input-output were antialiasing filtered by an 
eighth-order Bessel filter with a cut-off at 200 Hz and recorded at 1,000 Hz. Data was collected 
for a range of Mach numbers from M = 0.6 — 1.2 and dynamic pressures Q — 20 — 65 psf 
(pounds per square foot). 
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3.2 Data Analysis 


In this study, the objective was to develop a parsimonious model for control synthesis to improve 
gust-load alleviation and ride quality. As such, we focused our identification efforts on data 
collected at Mach 0.80, 0.95, 1.10 and Q = 30, 55, 60, 65 psf. At each Mach we used Q = 30 
psf as the estimation data and Q = 55, 60, 65 psf as the cross-validation data. At M = 0.95 
data was not collected at Q = 65 psf and, therefore, the estimated model could only be cross- 
validated at Q = 55 and 60 psf. This approach of modelling and validation was taken to assess 
the model’s predicative capability and robustness. 

Data was preprocessed to remove the linear trend, mean and outliers. The preprocessing step 
ensured that all unwanted low-frequency disturbances, offsets, trends and drifts were removed 
and allowed for an accurate representation of the system dynamics. 

3.3 Identification Procedures 

The identification process was performed in four stages to assess the full range of models dis- 
cussed in §2. In all cases the model order was estimated using the FPE, AIC and MDL ap- 
proaches. The optimal model order was deemed as one that produced the lowest prediction 
error of the aforementioned techniques. Once the model order was fixed, prediction error iden- 
tification (PEI) was used to estimate the unknown parameters and its predictive capability was 
computed for a cross-validation set as described previously (see §2.3). 

1. The ARX structure’s ability to model wind-tunnel data was assessed because it is the sim- 
plest input-output model and invoking the principle of parsimony. For the ARX model the 
order was determined by preselecting a range of model orders to search over. Specifically, 
n a =2-5, rib =1-5 and rik =1-10 was chosen as the search range, where n k denotes input 
delay. This range was selected to allow sufficient model complexity whilst maintaining 
an efficient system description. 

2. Further complexity was added to the ARX model to assess whether the same model struc- 
ture could account for significantly greater output variance with n a =2-10, rib =1-10 and 

rife =1-20. 

3. Using the best fit ARX model order to fix n a and rib we searched for an ARMAX order, 
namely, n c =1-20 that provided the smallest prediction error as the optimal ARMAX 
order. 
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4. The same procedure was followed to assess whether the B J model structure could account 
for more of the output variance than the ARMAX model with rid =1-20. 

4 RESULTS 
4.1 ARX Model 

Two ARX model orders were evaluated. Specifically, a model with order n a =5 , rib =1, ’>>k =1 
and n a =10 , rib =1, =1 was used for model development. The model order was estimated 

as discussed in §2.2. This analysis allowed for us to assess whether adding complexity to the 
same structure could account for significantly more of the output variance. The models were 
estimated using wind-tunnel data measured at Mach 1.10, Q — 30 psf which yielded structures 
of the form 


ARX5: 

A(q)y(n) 

= B(q)u(n) + e(n) where 

(13) 

A(q) 

— 1 H“ Qj\C[ i + &2Q ^ ^ ^ + cl§(J ° and 


B(q) 

= hq~ l 


ARX 10: 

A(q)y(n) 

= B(q)u(n ) + e{n) where 

(14) 

A(q) 

= 1 + di q 1 + ci2q 2 + (13 q 3 + CLiq 1 + ci^q J 
+ aQ q 6 + a^q 6 + a^q 6 + a^q 6 + dio 7 10 


B(q ) 

= hq- 1 . 



Figure 2 shows representative results for the ARX structure’s (Eqns.13 & 14) ability to represent 
S 4 T wind-tunnel data by evaluating it with cross-validation data at Mach 1.10, Q =55, 60, 65 
psf. This figure compares %QF of the predicted output for the ARX5 and ARX 10 models 
superimposed on top of measured data. The %QF’s obtained for the ARX5 model at 0 =55, 
60, 65 psf. are 82.97%, 84.94% and 84.94%, respectively. For the ARX10 model the %QF’s 
are 86.61%, 88.24% and 88.05%. 

Figure 2 demonstrates that adding complexity to the ARX model improves the %QF. As such, 
we deem that the ARX 10 model (Eqn. 14) is the better fit model due to its predicative capability. 


10 



Cross-Validation Fit: Mach 1.10, Q=55 


Measured 

ARX10 86.61 %QF 

ARX5 82.97%QF 


0.5 


O 0 



0 20 40 60 80 100 

Time (s) 


(a) 


0 


Cross-Validation Fit: Mach 1.10, Q=60 


Cross-Validation Fit: Mach 1.10, Q=65 



0 20 40 60 80 100 

Time (s) 


0 20 40 60 80 100 

Time (s) 


(b) 


(c) 


Figure 2: Measured and predicted output for ARX5 and ARX10 models at Mach 1.10. (a): Q = 55 psf. (b): 
Q = 60 psf. (c): Q = 65 psf. Solid line (“ — ”) measured output. Dash-dash line (“ — ”) predicted 
ARX10 model output. Dot-dot line (“...”) predicted ARX5 model output. 
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4.2 ARMAX Model 


Next, we used the ARX10 model oder as a starting point to develop a more complex ARMAX 
model to assess whether it could account for more output variance whilst maintaining an effi- 
cient model description. The model was estimated at Mach 1.10, Q = 30 psf and yielded a 
structure of the form 

ARMAX: 

A(q)y(n) = B(q)u(n) + C(q)e(n) where (15) 

C(q) = 1 + Ciq^ 1 + c 2 q~ 2 + c 3 q~ 3 + c A q~ 4 + c 5 g~ 5 
+ c 6 q~ 6 + c 7 q~ 7 + c 8 q~ 8 + c 9 q~ 9 + c 10 q~ 10 . 

Figure 3 shows representative results for the ARMAX structure’s (Eqn. 15) ability to represent 
S 4 T wind-tunnel data by evaluating it with cross-validation data at Mach 1.10, Q =55, 60, 65 
psf. The figure illustrates %QF of the predicted output for the ARMAX model superimposed 
on top of measured data. The %QF’s obtained for the ARMAX model at Q =55, 60, 65 psf. 
are 87.38%, 88.86% and 88.67%, respectively. Notice that although this model structure adds 
complexity, it is accounts for incrementally more of the output variance. 

4.3 Box- Jenkins Model 

Lastly, we used the ARMAX model (Eqn. 15) oder as a starting point to develop a more complex 
BJ model to assess whether it could account for more output variance whilst maintaining an 
efficient model description. The model was estimated at Mach 1.10, Q = 30 psf and yielded a 
structure of the form 


BJ: 



y(n) 

B{q) , \ 

= ^ — -u(n) + 

C(q ) 

— eyn) where 

F ( q ) 

D{q) 

m 

~ A(q) and 


D(q) 

— 1 + d\q 1 + 

d 2 q 2 + d 8 q 3 + d^q 1 + d§q 


+ d^q 6 + d 7 q 

7 + d 8 q 8 + dgq 9 + di^q 10 . 


Figure 4 shows representative results for the BJ structure’s (Eqn. 16) ability to represent S 4 T 
wind-tunnel data by evaluating it with cross-validation data at Mach 1.10, Q =55, 60, 65 psf. 
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Cross-Validation Fit: Mach 1.10, Q=55 




-1.5 


20 40 60 80 100 

Time (s) 


(a) 



Figure 3: Measured and predicted output for ARMAX model at Mach 1.10. (a): Q = 55 psf. (b): Q = 60 psf. (c): 
Q = 65 psf. Solid line (“ — ■”) measured output. Dot-dot line (“...”) predicted ARMAX model output. 
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Cross-Validation Fit: Mach 1.10, Q=55 


1.5 


Measured 

BJ 87.1 7%QF 
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Time (s) 


(a) 



(b) 


(c) 


Figure 4: Predicted output of the BJ model at Mach 1.10 superimposed on top of measured output, (a): Q = 55 
psf. (b): Q = 60 psf. (c): Q = 65 psf. Solid line (“ — ■”) measured output. Dot-dot line (“...”) predicted 
BJ model output. 
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The figure illustrates %QF of the predicted output for the BJ model superimposed on top of 
measured data. The %QF’s obtained for the BJ model at Q =55, 60, 65 psf. are 87.17%, 
88.70% and 88.48%, respectively. Although the BJ model offers greater complexity to model 
the observed data, it offers a slightly lower %QF. Using validation data, if the fit of a higher 
order deteriorates, it is an indication that the model complexity is too high [33]. 

Table 1 summarises the findings of our analysis. These results from Table 1 illustrate that model 



Dynamic 


Quality of Fit (%) 


Mach 

Pressure (psf) 

ARX5 

ARX 10 

ARMAX 

BJ 


55 

86.04 

87.63 

88.26 

88.05 

0.80 

60 

86.81 

88.25 

88.87 

88.61 


65 

86.16 

87.33 

87.79 

87.66 


55 

83.31 

85.56 

86.80 

86.05 

0.95 

60 

83.60 

85.97 

87.11 

86.36 


65 

X 

X 

X 

X 


55 

82.97 

86.61 

87.38 

87.17 

1.10 

60 

84.94 

88.24 

88.86 

88.70 


65 

84.94 

88.05 

88.67 

88.48 


Table 1: Summary of cross-validation results. Mach number and dynamic pressure versus model structure. 

fit increases with added complexity for the ARX model structure, improves incrementally for 
the ARMAX structure but decreases for the BJ. From these results we conclude that the higher- 
order ARX model may be sufficient for our purposes and as such we compare the predictive 
capability of this model to the FE based ASE model used for control design during wind-tunnel 
test [38,39]. 

4.4 Comparison of Tenth-Order ARX and ASE Models 

Lastly, we compared the ARX10 (Eqn. 14) and ASE model’s ability to predict measured data. 
The original FE based model had 17,196 degrees of freedom. Using a modal approach an 
eightieth-order ASE model was developed and used for comparison [38,39]. The ASE model 
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contained 60 structural and 20 unsteady aerodynamic lag states. 


Figure 5 illustrates results of how accurately, as %QF, the two models correspond to wind-tunnel 
data at Mach 0.80, 0.95, 1.1 and, Q =55, 60, 65 psf. The figure shows the predicted output for 


Cross-Validation Fit: Mach 0.80, Q=55 Cross-Validation Fit: Mach 0.80, Q=60 Cross-Validation Fit: Mach 0.80, Q=65 
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(e) 


Cross-Validation Fit: Mach 1.10, Q=55 Cross-Validation Fit: Mach 1.10, Q=60 Cross-Validation Fit: Mach 1.10, Q=65 
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Figure 5: Predicted outputs of the ASE and ARX10 model superimposed on top of measured output, (a-c): Mach 
0.80, Q = 55-65 psf. (d-e): Mach 0.95, Q = 55 - 60 psf. (c): Mach 1.10, Q = 55 - 65 psf. Solid 
line (“ — ”) measured output. Dot-dot line (“...”) predicted ARX10 model output. Dash-dot line 
predicted ASE model output. 


the ARX10 and ASE models superimposed on top of measured data. The %QF’s obtained for 
the ARX10 model, at all Mach numbers and Q’s, ranges from 85% to 88%. For the ASE 
model the %QF’s, for all Mach numbers and Q’s, ranges from 2% to 7%. Although the FE 
ASE has many more degrees of freedom, it was not able to provide better predictive capability 
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than the data-driven ARX model. Clearly, the tenth-order ARX model obtained using system 
identification methodology outperforms the FE ASE model. 


5 DISCUSSION 

This study explored the utility of system identification techniques to develop robust models 
for control synthesis. Initially, at each Mach number, an ARX model with maximum order 
of five was posed to the AIC, FPE and MDL techniques to estimate optimal lag and dynamic 
order. Next the ARX models were allowed greater flexibility with a maximum order of ten to 
assess whether the same structure could explain more of the output variance. Analysis of these 
results indicated that the tenth-order ARX models provided sufficient improvement in model fit 
to justify the added complexity. Therefore, at each Mach number the tenth-order ARX models 
were used as a starting point for ARMAX model development. The best fit ARMAX models 
were then used to develop BJ models. 

The results show that whilst the ARMAX models provide an improved fit, it was less than 1% 
for each case. The predictive capability of the BJ models was less than that of the ARMAX 
structures, which is a symptom of too much complexity. As such, these results indicate that 
the tenth-order ARX models were optimal in terms of parsimony and predictive capability to 
describe the dynamics of recoded wind-tunnel data. 

A comparison of the tenth-order ARX and 80th-order ASE model’s ability to explain the output 
variance revealed that the more complex ASE model was not able to achieve as high of a fit. The 
ARX model was able to attain a higher fit because the model was developed using measured 
data which often contains dynamics that are not captured in a FE based model. Conversely, the 
ASE model was developed from idealised assumptions about mass, damping, etc. which are 
idealisation and often do not closely hold with observations under experimental conditions. 

Although a different ASE model was developed for each Mach and Q, they yielded significantly 
reduced fits to measured data compared with the tenth-order ARX model. Notice that the ARX 
model at each Mach number was developed with estimation data at Q = 30 psf but was able to 
accurately predict the measured data at higher Q’s, namely, Q — 55 — 65 psf. This is the power 
of utilising a data-driven approach to model development. 

Often a model which yields a minimum system description yet provides good predictive ca- 
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pabilities is deemed as the best or optimal model. This is a trade-off between the ability to 
describe the system behaviour and parsimony. If two models can describe the system behaviour 
almost equally well, why choose the more complex one? For example, in -control synthesis 
the order of the controller is equal to the order of the model plus the order of the performance 
weights. Of course, it is not as computationally demanding to compute a controller of lower 
order. Moreover, selecting a model with greater complexity may lead to numerical issues if 
the model order is too high. Hence, minimum complexity often renders control synthesis more 
tractable [33,47]. 

A future study will utalise these data-driven models to develop control laws to assess model 
and controller design robustness. Initially, this study will be performed in a simulation environ- 
ment but with data collected from the wind-tunnel tests described in this paper. The successful 
demonstration of this work may lead to comparison of the robustness of control law design 
based on the two approaches, FE ASE and system identification models, on a supersonic flight 
test vehicle. 

6 CONCLUSION 

This study demonstrates the application of system identification techniques to develop parsi- 
monious and robust models directly from data with excellent predictive capability. The results 
show that a data-driven model is capable of predicting observed data for a larger operating 
point. This robust predictive power was demonstrated with models developed at Q = 30 psf 
yet able to accurately predict the measured output at higher Q’s. This superior predictive capa- 
bility allows for simpler control solutions and reduced modelling effort whilst traditional ASE 
based control strategies rely on computationally expensive models with little predictive power 
rendering control law design more expensive. 
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Introduction 


• Undesirable ASE interactions are a major concern in modern aircraft design 

- ASE interactions between aircraft structure, aerodynamics and flight control systems 
can lead to divergent oscillations resulting in catastrophic failure 

- Accurate models allow for robust control design - critical for aircraft safety, gust-load 
alleviation, ride quality, etc. 

• Finite element (FE) based models used in design process of aircraft to aid 
in the description of complex elastic and structural components 

- FE based models that accurately characterise the aerodynamic and structural compo- 
nents are of very high order (e.g. thousands of degrees of freedom) and computationally 
intensive 

• Fed to considerable activity in the areas of model and controller reduction 
techniques - reduced order model (ROM) techniques 
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ROM Techniques 




• Proper Orthogonal Decomposition (POD) (i.e. Karhunen-Loeve Procedure) 

- POD technique well known in statistical literature as principal-component analysis 

- Uses synthetic data from a high fidelity FE based model to capture the dominant char- 
acteristic information utilising an orthogonalisation process 

• Trajectory Piecewise Linearisation (TPWL) 

- Developed because standard POD procedures less robust for nonlinear problems 

- TPWL technique combines reduced-order modelling with linearization of the govern- 
ing equations as a solution to this problem 

• ROM techniques offer significant reduction of full FE based model, how- 
ever, often still too large for efficient control design 

• An alternative approach is data driven techniques to let the data dictate what 
the optimal model should be - i.e. system identification 


© 


IFASD 201 1 - Paris, France 


9 

Dryden Flight Research Center 


The Identification Problem 



Measured Input u(t) Measured Output z(t) 
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Identification of Stochastic or Noise Model 



Measured Input u(t) Measured Output z(t) 


- Relationship between w\{t) and z(t), given only system output, z(t) 

- u(t) assumed zero or constant 

- Commonly known as time-series analysis 

- Economic analysis, geophysical or astronomical phenomena, biological data 
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Identification of the Deterministic Model, Q 

^ ■ ■ • 



Measured Input u(t) Measured Output z(t) 


- Relationship between u(t) and y[t) - assumes w\(t) = 0 

- Input/output corrupted by noise, e u (t) & e z (t) - commonly assumes e u (t) = 0 

- Pursued when objective is to gain insight into the functioning of a system 

- Automotive industry, chemical plants, pulp & paper, biomedical modelling 
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Identification of Stochastic & Deterministic Models 



Measured Input u(t) Measured Output z(t) 


- Both input/output signals available for identification 

- Used when accurate predictions are desired 

- Design of model-based control systems for aircraft, spacecraft or robotics 
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Nonparametric Methods 


• Advantage 

- Provide convenient, robust means of characterising the dynamics of systems without 
requiring a priori assumptions regarding the system structure 


• Disadvantage 

- Nonparametric estimates of dynamics are difficult to relate to the structure and param- 
eters of the underlying physiological system 
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Parametric Methods 


• Disadvantage 

- Generally require a priori assumptions about the system order 


• Advantage 

- Provides a concise description of the system dynamics 

- Yield results that may be related directly to the system structure 
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Objective 


• Investigate two common approaches to model development for robust con- 
trol synthesis in the aerospace community 


• Reduced order aeroservoelastic modelling based on structural finite-element 
and computational fluid dynamics based aerodynamic models 


• Data-driven system identification procedure 
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Identification Procedure 


• Choose Model Structure 

- Determine model form and select parameters to include in a model 

• Model Order Selection 

- Determine number of input, output and error lags 

• Parameter Estimation 

- Determine values of unknown parameters 

• Model Validation 

- Assess whether identified nominal model can reproduce data from a plant 
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Model Structures 


. ARX Model 


y(n) 

y(n) 


-aiy{n - 1) a na y(n - 

+e(n) 

G(q)u(n) + H{q)e(n) where 


n a ) + biu(n — 



1) + b 2 u(n — 2) H h b nb u(n 


H{q) — 


1 

A {q) 



• ARMAX Model 



-aiy(n - 1) «n 0 2/(« - n a ) + b x u(n 


+e(n) + cie(n 

where G(q) = 


-1) + 
m 

a(qY 


+ c nc e(n - n c ) 

C(q) 


H(q) = 


Mq) 


1) + b 2 u(n — 2) H h b nb u(n - n b ) 


• Box-Jenkins Model 


y(n) = 

where 


m 

F(q) 

F(q) 


, - C(q) 

u(n) + W) 

= i + hq l 


e(n ) 

+ • • • + fn f Q Uf 


and D(q) 


1 + d\q 1 + • • • + d Ud q n<1 
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Model Order Selection 


Final Prediction Error (FPE) 


FPE = V 1 + 


2 p 


N — p 

where N is the number of data points, V is the prediction error, or the residual sum of 
squares and p is the number of model parameters 


Akaike’s Information Criterion (AIC) 


AIC = V 1 + 


2 p 

N 


Minimal Description Length (MDL) 


MDL = V 1 + 


pin N 
N 
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Cross-Validation 


ptt; 


• Model validation concerned with assessing whether a given nominal model 
can reproduce data from a plant, collected after some initial experiments to 
obtain estimation data 


• Model goodness is assessed by 



L _ w EjL (yn - y«) 2 \ 

V W E N nM 2 ) 


x 100 
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Experimental S 4 T Wind-Tunnel Data from NASA Langely Research Center 



• Data analysed used horizontal tail position input and structural accelerom- 
eter response output 
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Data Collection 


ptt; 


• Data recorded during transonic clearance 


• Aircraft model was perturbed with a log sine sweep input and had a mean 
value of 3.5 deg and ± 0.3 deg amplitude 


• Inputs applied to horizontal tail and output data collected from sensor lo- 
cated at nacelle inboard aft position 


• Data was collected for M = 0.6 — 1.2 and dynamic pressures Q = 20 — 65 
psf 
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Data Analysis 


• Objective was to develop a parsimonious model for control synthesis to 
improve gust-load alleviation and ride quality 


• Focused our identification efforts on data collected at Mach 0.80, 0.95, 1.10 
and Q = 30, 55, 60, 65 psf 


• Used Q = 30 psf as estimation data and Q = 55,60,65 psf as cross- 
validation data 


• At M = 0.95 data was not collected at Q = 65 psf and, therefore, the 
estimated model could only be cross-validated at Q = 55 and 60 psf. 
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Identification Procedures 


• Model order was estimated using FPE, AIC and MDL 

• Prediction error identification (PEI) used to estimate unknown parameters 

• Identification process was performed in four stages: 

1. ARX structure’s ability to model wind-tunnel data assessed because it is the simplest 
input-output model and invoking the principle of parsimony. Model order determined 
by preselecting a search range of: n a =2-5, rib =1-5, n k =1-10. Range selected to 
allow sufficient model complexity whilst maintaining an efficient system description. 

2. Complexity added to ARX model to assess whether the same model structure could 

account for significantly greater output variance with n a =2-10, rib =1-10, =1-20. 

3. Using best fit ARX model order to fix n a and rib we searched for an ARMAX order, 
namely, n c =1-20 that provided the smallest prediction error as the optimal order. 

4. The same procedure was followed to assess whether the BJ model structure could ac- 
count for more of the output variance than the ARMAX model with rid =1-20. 


© 


IFASD 201 1 - Paris, France 


9 

Dryden Flight Research Center 


Results: ARX Model.... 


• Two ARX model orders evaluated: n a =5, =1, =1 and n a =10, 

rib =1, rik =1 with estimation data from Mach 1.10, Q = 30 psf 


ARX5: 

/\ 

A(q)y(n) 

Mq) 

m 

ARX10: 


B(q)u(n) + e(n) where 

1 + diq + a 2 q^ 2 + <23 + a . and 

A 1 

biq 


A(q)y(n) = B(q)u(n) + e(n) where 

A(q) = 1 + diq" 1 + a 2 q~ 2 + d 3 t/~ 3 + d 4 (/~ 4 + a 5 <A 5 

+ a§q + a^q + a^q + a^q + a^q and 
B(q) = h q~ x . 
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....Results: ARX Model Cross-Validation at Mach 1.10, Q =55, 60, 65 psf 


Cross-Validation Fit: Mach 1.10, Q=55 


Measured 

ARX10 86.61 %QF 


ARX5 82.97%QF 
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Results: ARMAX Model.... 


• Used the ARX10 model order as a starting point to develop a more complex 
ARMAX model with estimation data from Machl.10, Q = 30 psf 


ARMAX: 

A{q)y{n) 

C(q) 


= B(q)u(n) + C(q)e(n) where 
= 1 + Cig 1 + c 2 q~ 2 + c 3 (A 3 + c 4 (A 4 + c 5 c/ 
+ c 6 q~ 6 + c 7 q~ 7 + c 8 q~ 8 + c 9 q~ 9 + Ci 0 q~ 10 
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....Results: ARMAX Model 


Cross-Validation Fit: Mach 1.10, Q=55 


1 .5 r 


- Measured 
ARMAX 87.38%QF 


C5 0 
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Cross-Validation Fit: Mach 1.10, Q=60 Cross-Validation Fit: Mach 1.10, Q=65 
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Results: Box- Jenkins Model.... 


• Used the ARM AX model order as a starting point to develop a more com- 
plex BJ model with estimation data from Machl.10, Q = 30 psf 


BJ: 

yip) 

f(q) 

D(q) 


B(q) C(q) 

= — — uin) + — — -ein) where 
F(q) D(q) 

~ A(q) and 

= 1 + d\q~ l + d 2 q~ 2 + d 3 g~ 3 + d^q~ A + d$q ~ 5 
+ d§q~ Q + d^q~ 7 + d%q~ 8 + d$q~ 9 + dioq~ 10 . 
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....Results: Box- Jenkins Model 


Cross-Validation Fit: Mach 1.10, Q=55 


1 .5 r 


C5 0 


-0.5 


-1.5 


- Measured 
BJ 87.17%QF 
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Cross-Validation Fit: Mach 1.10, Q=60 Cross-Validation Fit: Mach 1.10, Q=65 
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Summary of Cross-Validation: Mach & Dynamic Pressure Vs. Model 

— : > 


• Results illustrate that model fit increases with added complexity for the 
ARX model structure, improves incrementally for the ARMAX structure 
but decreases for the BJ 


• We conclude that the higher-order ARX model may be sufficient for our 
purposes and as such we compare the predictive capability of this model to 
the FE based ASE model used for control design during wind-tunnel test 
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Comparison of Tenth-Order ARX and ASE Models 


• Compared the ARX10 and ASE model’s ability to predict measured data 


• FE based model had 17,196 degrees of freedom, using a modal approach 
an eightieth-order ASE model was developed and used for comparison 


• ASE model contained 60 structural and 20 unsteady aerodynamic lag states 
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Predicted Outputs of ASE, ARX10 Models Measured Output 




Cross-Validation Fit: Mach 0.95, Q=55 
1 -5 Measured 


ARX1 0 85.56%QF 
ASE Model 1 ,638%QF 
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Cross-Validation Fit: Mach 1.10, Q=55 Cross-Validation Fit: Mach 1.10, Q=60 Cross-Validation Fit: Mach 1.10, Q=65 
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Conclusions 




• Study demonstrates the application of system identification techniques to 
develop parsimonious and robust models directly from data with excellent 
predictive capability 


• Results show that a data-driven model is capable of predicting observed 
data for a larger operating poin 


• Robust predictive power was demonstrated with models developed at Q = 
30 psf yet able to accurately predict the measured output at higher Q’s 


• Superior predictive capability allows for simpler control solutions and re- 
duced modelling effort whilst traditional ASE based control strategies rely 
on computationally expensive models with little predictive power rendering 
control law design more expensive 
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